A cluster heat bath method on a quasi-one dimensional Ising model 
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We have proposed a cluster heat bath(CHB) method of a quasi-one dimensional Ising model and 
demonstrated that it reproduces the exact results of the two dimensional Ising model with axially 
anisotropic exchange interactions. The point of the method is to select one of the equilibrium spin 
configurations of each of the chains which are subjected by effective fields of surrounding spins. 
Thus the chains are always in their equilibrium state during the simulation and the d dimensional 
model effectively becomes as a d — 1 dimensional model of fictitious spins with their freedom of 2 N 
and the simulation time is drastically reduced, where N is the number of the spins on the chain. 



Recently, progress of the computational ability enable 
us to execute large scale computer simulations and more 
accurate nature of systems have been revealed. |ij How- 
ever, low temperature phases of complex systems such as 
spin glasses have not well recognized yet, because they 
have many local minima and huge CPU time is neces- 
sary to perform computer simulations. To remove the 
difficulty, various algorithms have been proposed such as 
cluster algorithms and extended ensemble methods. 

Similar problems also occur in studying the phase 
transition even in non- complex systems whose ground 
state is known. Among them, quasi-one dimensional 
magnetic systems are of great interest, because many at- 
tractive substances have been found such as CSC0CI3, |t]] 
CsCoBr 3 H and CSMCI3 for which successive phase 
transitions are observed. However, so long as we know, 
no effective algorithm has been proposed yet to treat 
those systems. 

In this letter, we discuss the difficulty faced on Monte 
Carlo(MC) simulations of qausi-one dimensional systems 
and propose a new effective method to remove it. Our 
method is a cluster algorithm but different from previous 
ones [^],[| in which all the spins in the cluster are flipped 
simultaneously in accordance with the detailed valance 
condition. In our algorithm, in contrast with this, one of 
the possible spin configurations of the cluster is selected 
according with its Boltzmann weight. So we may call it 
as a cluster heat bath(CHB) method. In the qausi-one di- 
mensional system, this can be readily done, because the 
clusters may by chosen as one-dimensional chains and 
a recursion relation may be used to determine the spin 
configuration. 

We start with a two dimensional ferromagnetic Ising 
model with axially anisotropic nearest neighbor interac- 
tions described by the Hamiltonian: 
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where J and J'(<C J) are coupling constants along the 
longitudinal and the transverse directions, respectively, 
and Sg = ±1. This model was already solved by On- 
sager pfl] about 50 years ago. However, the MC simula- 
tion of this model is not so easy. Suppose that we make 
a MC simulation of the model on an M x N lattice. As 



the temperature is decreased from a high temperature, 
the spin correlation in every chain along the longitudinal 
direction develops first. The correlation length £(T) of 
the chain may be given as 



£(T) ~ exp(2 J/T), 



(2) 



where T is the temperature measured in a unit of ks = 1. 
As the temperature is decreased further, the correlation 
along the transverse direction develops and the long- 
range order will take place at T — T c in the thermo- 
dynamic limit. However, if N < £(T C ), the spins on 
each of the chains will order ferromagnetically even at 
T > T c and they will behave like a single-spin. So that 
the system will not show true critical behavior of the two 
dimensional model. Therefore, we must treat the lattice 
with N 3> £,{T C ) 3> 1 and much CPU time is necessary to 
make the MC simulation. If, however, there exists some 
algorithm which leads an equilibrium spin configuration 
of individual chain, we may reduce the CPU time, be- 
cause a large MC step is wasted to equilibrate the spin 
configuration of every chain. Here we propose that algo- 
rithm. 

We consider the two dimensional Ising model of eq.(l) 
on a rectangle lattice of M x N. For simplicity, we im- 
pose the periodic boundary condition to the transverse 
direction, and the open boundary condition to the longi- 
tudinal direction. We consider the chain of N spins along 
the longitudinal direction whose Hamiltonian is given by 
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hi = J' \Sij+i + Sij—i), 
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where <7j = Sij. Then the problem is to determine the 
spin configuration of the one dimensional chain with ran- 
dom fields hi. This can be readily done as follows. 
We define the function: 

fn(o-) = ^ exp{if (CT1CT2 + 0-20-3 H h o„_io) 

<Tl f O"2,-",0 r n _l 

+ (Cioi + CW 2 + • • • + C n a)}, (5) 

where K = J/T and C n = h n /T. f n (a) can be readily 
obtained from the recursion formula 
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/„_i(CT„_i)exp{i\cr„_ier + C n o} n>2, 
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with fi(a) — exp(Cicr). Then the probabilities Pn(&n — 
±1) of the end spin being on = ±1 are given as 



P N (a N = ±1) = 



fN(cr N = ±1) 



(7) 



where Zn(= /jv(1) + /n(~ 1)) is the partition function 
of the chain. Then we can determine the value of on by 
using a uniform random number between to 1. The 
next step is to determine the value of on-i under the 
condition that the value of on is fixed either on = 1 or 
(7 j\[ = — 1 . We can readily do it by changing the effective 
field on <jn-i from Kn—x to hN-i + JoN and the function 
/tv-i(o-) to 

/jV-l(o-) = ^ fN-2((JN-2) 

<XjV-2 

X exp{i^(TAr_2C r + Cn-IC + KoOn} 

= fN-i(<^)exp(Kaa N ). (8) 
Then the probabilities Pn—i(&n-i = il) are given as 



-Pjv-i(o-at-i = ±1) 



/jv-i(ctv-i 



±1) 



(9) 
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where 



^N-x — ^ /jv-i(cat-i) 

= /Ar(crjv)exp(-CjvCTjv), 



(10) 



and the value of on-i can be determined. Repeating 
this procedure from iV - 1 to 1, we can determine the 
spin configuration of the chain. Using this algorithm, the 
spin configuration of the lattice is updated from chain to 
chain. 

Note that the spin configuration of the chain 
f (ci,<72, ...,ctjv) realized in this procedure accords with 
the Boltzmann weight. That is, 
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So it is natural to call this method as a chain heat bath 
method or more generally a cluster heat bath method and 
to abbreviate it to a CHB method. 
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FIG. 1. The relaxation processes in the lattice with 50 x 500 
at (a)T = 0.6 J and (b)T = 0.4J. The broken line is the exact 
value by Onsager(ref.lO). 



We apply the CHB method to the model (1) with 
J' = 0.01J. The transition temperature is obtained from 
eq.(17) in ref.10 as T c ~ 0.51J. We examine the relax- 
ation of the energy starting from a random spin configu- 
ration. In Figs. [j](a) and (b), we present typical results 
of the MC step dependence of the energy at temperatures 
above and below T c , respectively, together with those by 
a conventional heat bath(HB) method. We find that the 
CHB method reproduces the exact result within a MC 
step much smaller than that of the HB method. We 
need only 1 or 2 MC steps to get the equilibrium value 
for T > T c and several dozens of MC steps for T < T c , 
whereas we could not get it within 10000 MC steps in the 
HB method even for T > T c . In Figs. |(a) and (b), we 
present snapshots of the spin configuration for T < T c at 
different MC steps, t, in the HB and the CHB methods, 
respectively. In the CHB method, even at t = 2, every 
chain seems to reach its equilibrium state, and domains 
are formed. Smaller domains shrink and disappear at 
t ~ 35. On the other hand, in the HB method, even 
at t = 1000, many of the chains are still not in their 
equilibrium state. These facts clearly reveal that the fast 
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relaxation in the CHB method comes from the instanta- 
neous relaxation of every chain. 

(a) HB 
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(b) CHB 




2 step 8 step 32 step 

FIG. 2. Snapshots of the spin configuration in the lattice 
with 50 x 500 at T = 0.4 J, where the longitudinal lattice space 
is reduced to fifty percent. Up and down spins are described 
by black and gray points, respectively. 

We tentatively define the relaxation time r as the step 
at which the energy firstly reaches the exact value. We 
made five runs starting from different random spin con- 



figurations and estimated it. In Fig. ||, we plot their 
average (r) and standard deviation At as functions of 
the temperature. 
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FIG. 3. Temperature dependences of the relaxation time 
(t). The symbols of square and diamond mean those for lat- 
tices with 26 x 250 and 50 x 500, and those of filled and open 
ones mean those in the HB and the CHB methods, respec- 
tively. 

In the HB method, (r) increases rapidly as the temper- 
ature is lowered towards T c and a large size dependence 
is seen. In the CHB method, for T > T c , (r) depends lit- 
tle on both the temperature and the lattice size. It once 
increases at T ~ T c then decreases as the temperature is 
decreased further. Even at T ~ T c , it is less than several 
dozens of MC steps. 
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FIG. 4. The specific heat of the Ising model with 
J' = 0.01J on the square lattice. Open and filled symbols 
mean those in the lattices with 50 x 500 and 100 x 1000, respec- 
tively, and a dotted line is the exact result by Onsager(ref.lO). 

The specific heat calculated by the CHB method is 
shown in the Fig. |]. The result is good agreement with 
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the exact value by Onsager in the whole temperature 
range. Especially, the sharp peak at T = T c is repro- 
duced well. 

We have proposed a cluster heat bath(CHB) method 
of a quasi-one dimensional Ising model and demonstrated 
that it reproduces exact results of the two dimensional 
Ising model within a MC step much smaller than that of 
the conventional MC method. The point of the method 
is to select one of the equilibrium spin configurations 
in each of the chains which are subjected by effective 
fields of surrounding spins. That is, each of the chains 
is always in its equilibrium state during the simulation. 
Then the d dimensional model becomes a d — 1 dimen- 
sional model of fictitious spins with the freedom of 2 N 
and the simulation time is drastically reduced. In fact, 
we could simulate a quasi-one dimensional Ising model 
on 36 x 34 x 1200 hexagonal lattice which is a model 
of CSC0CI3 and CsCoBr3 within a reasonable computer 
time [pd| and explain a change in the spin structure with- 
out any anomaly of the specific heat observed at a low 
temperature |8| as well as a usual phase transition at a 
higher temperature. 

Finally, we should note that we could readily apply the 
CHB method to usual Ising models by treating columns 
of the Ising spins of K x L x N(K, L = 1, 2, 3, 4, ...) us- 
ing a transfer matrix technique. The CHB method was 
found to be also effective for complex Ising models such 
as Ising spin glasses. Results will be reported elsewhere 
in a separate paper. Jllfl 
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